rm(list=ls())
library(vegan)
#library(SpiecEasi)
library(psych)##
library(igraph)
#setwd("~/GaoC/R")
library(beepr)
setwd("/Users/chenggao/Google_Drive/EPICON.NC.2022")
load("Bac.Fung.data.prep.stopby8datasets.amfL.RSLZxBF.Rarefaction.2021.08.13.rdata")
BF0<-read.csv("Microbiome/0000BacteriaFungi.1029x1293.13165.2019.07.10.csv", head = T, row.names = 1)
ID.tmp<-BF0[,c("Kingdom", "Kingdom1", "Phylum", "Class", "Order", "Family", "Genus", "Funguild", "OTU.ID", "Morph", "Morph1")]
#habitat <- "Root"; treatment <- "Control"; tpA <- 2; tpB <- 9
fq <- 7; abu <- 19;
SpMan<-function(fq, abu, treatment, tpA, tpB){
d1 <- d1.raw [env.R$Treatment1==treatment & env.R$TP > tpA & env.R$TP < tpB , ]
d1 <- d1 [,specnumber(t(d1)) > fq & colSums(d1)> abu]
d2 <- d2.raw [env.R$Treatment1==treatment & env.R$TP > tpA & env.R$TP < tpB, ]
d2 <- d2[, specnumber(t(d2)) > fq & colSums(d2)> abu]
env.R.tmp0 <- env.R [env.R$Treatment1==treatment & env.R$TP > tpA & env.R$TP < tpB , ]
d3 <- d3.raw [env.Z$Treatment1==treatment & env.Z$TP > tpA & env.Z$TP < tpB , ]
d3 <- d3 [,specnumber(t(d3)) > fq & colSums(d3)> abu]
d4 <- d4.raw [env.Z$Treatment1==treatment & env.Z$TP > tpA & env.Z$TP < tpB , ]
d4 <- d4 [,specnumber(t(d4)) > fq & colSums(d4)> abu]
env.Z.tmp0 <- env.Z [env.Z$Treatment1==treatment & env.Z$TP > tpA & env.Z$TP < tpB , ]
d5 <- d5.raw [env.S$Treatment1==treatment & env.S$TP > tpA & env.S$TP < tpB, ]
d5 <- d5[, specnumber(t(d5)) > fq & colSums(d5)> abu]
d6 <- d6.raw [env.S$Treatment1==treatment & env.S$TP > tpA & env.S$TP < tpB , ]
d6 <- d6 [,specnumber(t(d6)) > fq & colSums(d6)> abu]
env.S.tmp0 <- env.S [env.S$Treatment1==treatment & env.S$TP > tpA & env.S$TP < tpB , ]
d7 <- d7.raw [env.L$Treatment1==treatment & env.L$TP > tpA & env.L$TP < tpB , ]
d7 <- d7[,specnumber(t(d7)) > fq & colSums(d7)> abu]
d8 <- d8.raw [env.L$Treatment1==treatment & env.L$TP > tpA & env.L$TP < tpB, ]
d8 <- d8[, specnumber(t(d8)) > fq & colSums(d8)> abu]
env.L.tmp0 <- env.L [env.L$Treatment1==treatment & env.L$TP > tpA & env.L$TP < tpB , ]
saveRDS(env.L.tmp0, paste0("Fig.2/BFT.SpMan.df.2022.05.31/data.env.Leaf.", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
saveRDS(env.S.tmp0, paste0("Fig.2/BFT.SpMan.df.2022.05.31/data.env.Soil.", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
saveRDS(env.Z.tmp0, paste0("Fig.2/BFT.SpMan.df.2022.05.31/data.env.Rhizosphere.", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
saveRDS(env.R.tmp0, paste0("Fig.2/BFT.SpMan.df.2022.05.31/data.env.Root.", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
saveRDS(d1, paste0("Fig.2/BFT.SpMan.df.2022.05.31/data.bac.Root.", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
saveRDS(d2, paste0("Fig.2/BFT.SpMan.df.2022.05.31/data.fung.Root.", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
saveRDS(d3, paste0("Fig.2/BFT.SpMan.df.2022.05.31/data.bac.Rhizosphere.", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
saveRDS(d4, paste0("Fig.2/BFT.SpMan.df.2022.05.31/data.fung.Rhizosphere.", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
saveRDS(d5, paste0("Fig.2/BFT.SpMan.df.2022.05.31/data.bac.Soil.", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
saveRDS(d6, paste0("Fig.2/BFT.SpMan.df.2022.05.31/data.fung.Soil.", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
saveRDS(d7, paste0("Fig.2/BFT.SpMan.df.2022.05.31/data.bac.Leaf.", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
saveRDS(d8, paste0("Fig.2/BFT.SpMan.df.2022.05.31/data.fung.Leaf.", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
d12<-cbind(d1,d2)
spman.d12 = corr.test(d12, use="pairwise",method="spearman",adjust="fdr", alpha=.05, ci=FALSE)
saveRDS(spman.d12, paste0("Fig.2/BFT.SpMan.df.2022.05.31/SpMan.crossBF.Root.", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
d34<-cbind(d3,d4)
spman.d34 = corr.test(d34, use="pairwise",method="spearman",adjust="fdr", alpha=.05, ci=FALSE)
saveRDS(spman.d34, paste0("Fig.2/BFT.SpMan.df.2022.05.31/SpMan.crossBF.Rhizosphere.", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
d56<-cbind(d5,d6)
spman.d56 = corr.test(d56, use="pairwise",method="spearman",adjust="fdr", alpha=.05, ci=FALSE)
saveRDS(spman.d56, paste0("Fig.2/BFT.SpMan.df.2022.05.31/SpMan.crossBF.Soil.", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
d78<-cbind(d7,d8)
spman.d78 = corr.test(d78, use="pairwise",method="spearman",adjust="fdr", alpha=.05, ci=FALSE)
saveRDS(spman.d78, paste0("Fig.2/BFT.SpMan.df.2022.05.31/SpMan.crossBF.Leaf.", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
beep()
}
#SpMan(fq, abu, "Control", 2, 9)
#SpMan(fq, abu, "Pre_flowering_drought", 2, 9)
#SpMan(fq, abu, "Control", 8, 18)
#SpMan(fq, abu, "Pre_flowering_drought", 8, 18)
#SpMan(fq, abu, "Control", 9, 18)
#SpMan(fq, abu, "Post_flowering_drought", 9, 18)
######################################
# all correlations #
#####################################
rm(list=ls())
fq <- 7; abu <- 19;
#setwd("~/GaoC/R")
setwd("/Users/chenggao/Google_Drive/EPICON.NC.2022")
BF0<-read.csv("Microbiome/0000BacteriaFungi.1029x1293.13165.2019.07.10.csv", head = T, row.names = 1)
ID.tmp<-BF0[,c("Kingdom", "Kingdom1", "Phylum", "Class", "Order", "Family", "Genus", "Funguild", "OTU.ID", "Morph", "Morph1")]
#habitat <- "Root"; treatment <- "Control"; tpA <- 2; tpB <- 9
#r.cutoff = 0.6
#p.cutoff = 0.05
SpMan.Rall<-function(fq, abu, habitat, treatment, tpA, tpB){
d1 <- readRDS(paste0("Fig.2/BFT.SpMan.df.2022.05.31/data.bac.",habitat,".", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
d2 <- readRDS(paste0("Fig.2/BFT.SpMan.df.2022.05.31/data.fung.",habitat,".", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
d12 <- data.frame(d1, d2)
spman.r0 <- readRDS(paste0("Fig.2/BFT.SpMan.df.2022.05.31/SpMan.crossBF.",habitat,".", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
Cor<-as.matrix(spman.r0$r)
Cor.df<-data.frame(row=rownames(Cor)[row(Cor)[upper.tri(Cor)]],
col=colnames(Cor)[col(Cor)[upper.tri(Cor)]], Cor=Cor[upper.tri(Cor)])
P0<-as.matrix(spman.r0$p)
P.df<-data.frame(row=rownames(P0)[row(P0)[upper.tri(P0)]],
col=colnames(P0)[col(P0)[upper.tri(P0)]], p=P0[upper.tri(P0)])
da.tmp<-df.sig<- df <- data.frame(Cor.df, Habitat = habitat, Treatment = treatment, TPA = tpA)
#df[ df$Cor > r.cutoff & df$p < p.cutoff,]
V1<-data.frame("v1"=da.tmp$row); V2<-data.frame("v2"=da.tmp$col)
IDsub1<-ID.tmp[ID.tmp$OTU.ID %in% V1$v1, ]; IDsub2<-ID.tmp[ID.tmp$OTU.ID %in% V2$v2, ]
V1$id <- 1:nrow(V1); V2$id <- 1:nrow(V2)
M1<-merge(V1, IDsub1, by.x = "v1", by.y = "OTU.ID", all.x= T); M1<-M1[order(M1$id), ]
M2<-merge(V2, IDsub2, by.x = "v2", by.y = "OTU.ID", all.x = T); M2<-M2[order(M2$id), ]
df.tmp<-data.frame(da.tmp, M1, M2)
df.tmp$Link <- NA
df.tmp$Link[df.tmp$Kingdom=="Eukaryote" & df.tmp$Kingdom.1=="Eukaryote"] <- "Fun_Fun"
df.tmp$Link[df.tmp$Kingdom != df.tmp$Kingdom.1] <- "Fun_Bac"
df.tmp$Link[df.tmp$Kingdom=="Prokaryote" & df.tmp$Kingdom.1=="Prokaryote"] <- "Bac_Bac"
saveRDS(df.tmp, paste0("Fig.2/BFT.SpMan.all.df.2022.05.31/SpMan.Rall.Bac-Fung.",habitat,".", treatment,".",tpA,".",tpB,".",fq,".",abu,".RDS"))
}
#SpMan.Rall(fq, abu,"Root", "Control", 2, 9)
#SpMan.Rall(fq, abu,"Root", "Pre_flowering_drought", 2, 9)
#SpMan.Rall(fq, abu,"Root", "Control", 8, 18)
#SpMan.Rall(fq, abu,"Root", "Pre_flowering_drought", 8, 18)
#SpMan.Rall(fq, abu,"Rhizosphere", "Control", 2, 9)
#SpMan.Rall(fq, abu,"Rhizosphere", "Pre_flowering_drought", 2, 9)
#SpMan.Rall(fq, abu,"Rhizosphere", "Control", 8, 18)
#SpMan.Rall(fq, abu,"Rhizosphere", "Pre_flowering_drought", 8, 18)
#SpMan.Rall(fq, abu,"Soil", "Control", 2, 9)
#SpMan.Rall(fq, abu,"Soil", "Pre_flowering_drought", 2, 9)
#SpMan.Rall(fq, abu,"Soil", "Control", 8, 18)
#SpMan.Rall(fq, abu,"Soil", "Pre_flowering_drought", 8, 18)
#SpMan.Rall(fq, abu,"Leaf", "Control", 2, 9)
#SpMan.Rall(fq, abu,"Leaf", "Pre_flowering_drought", 2, 9)
#SpMan.Rall(fq, abu,"Leaf", "Control", 8, 18)
#SpMan.Rall(fq, abu,"Leaf", "Pre_flowering_drought", 8, 18)
rm(list=ls())
library(tidyverse)
library(vegan)##decostand
#library(WGCNA)##corAndPvalue
library(multtest)##mt.rawp2adjp
library(igraph)##graph.adjacency delete.vertices
library(psych)##
library(brainGraph)##
library(ggrepel)
library(dplyr)
library(plyr)
library(readr)
#library(SpiecEasi) #clr
library(tidyverse)
library(ggthemes)
library(beepr)
setwd("/Users/chenggao/Google_Drive/EPICON.NC.2022")
mydir = "/Users/chenggao/Google_Drive/EPICON.NC.2022/Fig.2/BFT.SpMan.all.df.2022.05.31"
da.Rdis <- list.files(path=mydir, full.names=TRUE) %>% map_dfr(readRDS)
da.Rdis$TPA<-factor(da.Rdis$TPA, labels=c("Drought_period" , "Rewetting_period"))
da.Rdis$Treatment<-factor(da.Rdis$Treatment, labels = c("Control", "Drought/Rewetting"))
da.Rdis$Links<-factor(da.Rdis$Link)
da.Rdis$intera<-factor(paste(da.Rdis$Links, da.Rdis$Treatment, da.Rdis$TPA, da.Rdis$Habitat))
#levels(da.Rdis$intera)
da.Rdis <-droplevels(da.Rdis[!(da.Rdis$Habitat=="Leaf" & da.Rdis$TPA=="Drought_period" & da.Rdis$Links!= "Bac_Bac"), ])
#da.Rdis$Links
library(splitstackshape) #cSplit
dia_me <- ddply(da.Rdis, .(intera), numcolwise(median))
dia_med<-cSplit(dia_me, "intera", " ")
dia_med$Links<-dia_med$intera_1
dia_med$Treatment<-dia_med$intera_2
dia_med$TPA<-dia_med$intera_3
dia_med$Habitat<-dia_med$intera_4
library(forcats) # fct_collapse
#da.Rdis$Treatment2<-fct_collapse(da.Rdis$Treatment, Stress = c('PRE', 'POST'))
#dia_med$Treatment2<-fct_collapse(dia_med$Treatment, Stress = c('PRE', 'POST'))
da.Rdis$Habitat<-factor(da.Rdis$Habitat, levels = c("Root", "Rhizosphere", "Soil", "Leaf"))
dia_med$Habitat<-factor(dia_med$Habitat, levels = c("Root", "Rhizosphere", "Soil", "Leaf"))
#sum(is.na(da.Rdis$Habitat))
#colnames(da.Rdis)
## figure size 14 x 4, density plot
ggplot(da.Rdis, aes(x= Cor, color = Links, linetype = Treatment)) +
geom_density(size=0.5)+
scale_color_manual(values = c("black", "blue","red", "cyan", "yellow", "red"))+
scale_linetype_manual(values = c(3, 1))+
xlab("Spearman's Rho")+ylab("Density")+
guides(fill=guide_legend(title= "Type"))+
geom_vline(data=dia_med, aes(xintercept=Cor, colour=Links, linetype = Treatment), size=0.3) +
facet_grid(Habitat~TPA)+theme_tufte()+
theme(strip.text = element_text(size = 12,face="bold"),
panel.spacing = unit(0, "lines"),
legend.title = element_text(colour="black", size=12, face="bold"),
legend.text = element_text(colour="black", size=12, face="bold"),
axis.text=element_text(size=8,face="bold"),
axis.title=element_text(size=12,face="bold"))

#ggsave("BFT.bac.fung.RSLZ.Rdensity.2021.08.27.pdf", width = 16, height = 10)
beep()
ggplot(da.Rdis, aes(x= as.numeric(Treatment),y= Cor, color = Habitat, linetype= Links)) +
geom_smooth(method='lm', size =0.5, se = FALSE)+
scale_color_manual(values = c("black", "blue","red", "cyan","brown", "yellow"))+
scale_linetype_manual(values = c(3, 1, 2))+
xlab("")+ylab("All Spearman's Rho")+theme_bw()+
guides(fill=guide_legend(title= "Type"))+
facet_wrap(~TPA, ncol =2, strip.position = "top")+
scale_x_continuous(breaks = c(1, 2),labels = c("Control", "Stress"))+
theme(strip.text = element_text(size = 10,face="bold"),
panel.spacing = unit(0, "lines"),
legend.title = element_text(colour="black", size=10, face="bold"),
legend.text = element_text(colour="black", size=10, face="bold"),
axis.text.x=element_text(colour="black",size=10,face="bold",angle = 90, hjust = 1),
axis.text.y=element_text(colour="black",size=8,face="bold"),
axis.title=element_text(colour="black",size=10,face="bold"))

#ggsave("BFT.bac.fung.RSLZ.Rall.CK.Stress.2021.08.27.pdf", width = 6, height = 3)
library(beepr)
beep()
ggplot(da.Rdis[da.Rdis$Cor>0,], aes(x= as.numeric(Treatment),y= Cor, color = Habitat, linetype= Links)) +
geom_smooth(method='lm', size =0.5, se = FALSE)+
scale_color_manual(values = c("black", "blue","red", "cyan","brown", "yellow"))+
scale_linetype_manual(values = c(3, 1, 2))+
xlab("")+ylab("Postive Spearman's Rho")+theme_bw()+
guides(fill=guide_legend(title= "Type"))+
facet_wrap(~TPA, ncol =2, strip.position = "top")+
scale_x_continuous(breaks = c(1, 2),labels = c("Control", "Stress"))+
theme(strip.text = element_text(size = 10,face="bold"),
panel.spacing = unit(0, "lines"),
legend.title = element_text(colour="black", size=10, face="bold"),
legend.text = element_text(colour="black", size=10, face="bold"),
axis.text.x=element_text(colour="black",size=10,face="bold",angle = 90, hjust = 1),
axis.text.y=element_text(colour="black",size=8,face="bold"),
axis.title=element_text(colour="black",size=10,face="bold"))

#ggsave("BFT.bac.fung.RSLZ.Rpos.CK.Stress.2021.08.27.pdf", width = 6, height = 3)
library(beepr)
beep()
ggplot(da.Rdis[da.Rdis$Cor<0,], aes(x= as.numeric(Treatment),y= Cor * -1, color = Habitat, linetype= Links)) +
geom_smooth(method='lm', size =0.5, se = FALSE)+
scale_color_manual(values = c("black", "blue","red", "cyan","brown", "yellow"))+
scale_linetype_manual(values = c(3, 1, 2))+
xlab("")+ylab("Negative Spearman's Rho * -1")+theme_bw()+
guides(fill=guide_legend(title= "Type"))+
facet_wrap(~TPA, ncol =2, strip.position = "top")+
scale_x_continuous(breaks = c(1, 2),labels = c("Control", "Stress"))+
theme(strip.text = element_text(size = 10,face="bold"),
panel.spacing = unit(0, "lines"),
legend.title = element_text(colour="black", size=10, face="bold"),
legend.text = element_text(colour="black", size=10, face="bold"),
axis.text.x=element_text(colour="black",size=10,face="bold",angle = 90, hjust = 1),
axis.text.y=element_text(colour="black",size=8,face="bold"),
axis.title=element_text(colour="black",size=10,face="bold"))

#ggsave("BFT.bac.fung.RSLZ.Rneg.CK.Stress.2021.08.27.pdf", width = 6, height = 3)
library(beepr)
beep()
ggplot(da.Rdis, aes(x= as.numeric(Treatment),y= Cor, color = Habitat)) +
geom_smooth(method='lm', size =5, se = FALSE)+
scale_color_manual(values = c("black", "blue","red", "cyan","brown", "yellow"))+
#scale_linetype_manual(values = c(3, 1, 2))+
geom_jitter(size = 0.1, alpha = 0.1)+
xlab("")+ylab("All Spearman's Rho")+theme_bw()+
guides(fill=guide_legend(title= "Type"))+
facet_grid(Habitat~TPA+Links)+
scale_x_continuous(breaks = c(1, 2),labels = c("Control", "Stress"))+
theme(strip.text = element_text(size = 10,face="bold"),
panel.spacing = unit(0, "lines"),
legend.title = element_text(colour="black", size=10, face="bold"),
legend.text = element_text(colour="black", size=10, face="bold"),
axis.text.x=element_text(colour="black",size=10,face="bold",angle = 90, hjust = 1),
axis.text.y=element_text(colour="black",size=8,face="bold"),
axis.title=element_text(colour="black",size=10,face="bold"))

library(beepr)
beep()
ggplot(da.Rdis[da.Rdis$Cor>0,], aes(x= as.numeric(Treatment),y= Cor, color = Habitat)) +
geom_smooth(method='lm', size =5, se = FALSE)+
scale_color_manual(values = c("black", "blue","red", "cyan","brown", "yellow"))+
#scale_linetype_manual(values = c(3, 1, 2))+
geom_jitter(size = 0.1, alpha = 0.1)+
xlab("")+ylab("Postive Spearman's Rho")+theme_bw()+
guides(fill=guide_legend(title= "Type"))+
facet_grid(Habitat~TPA+Links)+
scale_x_continuous(breaks = c(1, 2),labels = c("Control", "Stress"))+
theme(strip.text = element_text(size = 10,face="bold"),
panel.spacing = unit(0, "lines"),
legend.title = element_text(colour="black", size=10, face="bold"),
legend.text = element_text(colour="black", size=10, face="bold"),
axis.text.x=element_text(colour="black",size=10,face="bold",angle = 90, hjust = 1),
axis.text.y=element_text(colour="black",size=8,face="bold"),
axis.title=element_text(colour="black",size=10,face="bold"))

beep()
ggplot(da.Rdis[da.Rdis$Cor<0,], aes(x= as.numeric(Treatment),y= Cor* -1, color = Habitat)) +
geom_smooth(method='lm', size =5, se = FALSE)+
scale_color_manual(values = c("black", "blue","red", "cyan","brown", "yellow"))+
#scale_linetype_manual(values = c(3, 1, 2))+
geom_jitter(size = 0.1, alpha = 0.1)+
xlab("")+ylab("Negative Spearman's Rho * -1")+theme_bw()+
guides(fill=guide_legend(title= "Type"))+
facet_grid(Habitat~TPA+Links)+
scale_x_continuous(breaks = c(1, 2),labels = c("Control", "Stress"))+
theme(strip.text = element_text(size = 10,face="bold"),
panel.spacing = unit(0, "lines"),
legend.title = element_text(colour="black", size=10, face="bold"),
legend.text = element_text(colour="black", size=10, face="bold"),
axis.text.x=element_text(colour="black",size=10,face="bold",angle = 90, hjust = 1),
axis.text.y=element_text(colour="black",size=8,face="bold"),
axis.title=element_text(colour="black",size=10,face="bold"))

beep()